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We develop a boundary element method to calculate Van der Waals interactions for systems 
composed of domains of spatially constant dielectric response. We achieve this by rewriting the 
interaction energy expression exclusively in terms of surface integrals of surface operators. We 
validate this approach in the Lifshitz case and give numerical results for the interaction of two 
spheres as well as the van der Waals self-interaction of a uniaxial ellipsoid. Our method is simple 
to implement and is particularly suitable for a full, non-perturbative numerical evaluation of non- 
retarded van der Waals interactions between objects of a completely general shape. 



The problem of van der Waals interactions between 
objects of a general, low symmetry, shape is difficult [l|,0| 
and is related to the problem of shape dependence of 
eigenvalues of the wave equation in finite domains that 
has been consistently formulated and solved only for very 
restrictive conditions [3j. For general geometries the van 
der Waals interaction energy has been obtained mostly 
either in terms of a perturbation expansion in geometric 
deviations from the case of high symmetry (see Q and 
references therein) or a perturbation expansion in the 
dielectric contrast (see [f| and references therein). 

The result of Golestanian [j| is particularly germane 
for our point of departure since it deals with van der 
Waals interactions in general geometries, ft is based on 
a path integral field formulation and expresses the van 
der Waals interaction as a perturbation series in the spa- 
tial contrast of the polarizability profile e — 1. A different 
approach, based on the correlation effect of the density 
functional theory, is proposed in [(| and leads to an ap- 
proximate result for van der Waals interaction for general 
geometries in terms of an expansion in terms of 1 — e _1 . 
Both these perturbation results lead to tracelog formulas 
that contain volume integrals across the whole space. 

Below we will give an easily implementable numerical 
method for calculating non-retarded Van der Waals inter- 
actions for systems that are comprised of general-shape 
domains of spatially constant dielectric response. The 
expression we derive here is applicable to a wide class of 
geometries and is based on a surface trace reformulation 
of the interaction energy 0] . It is not based on any series 
expansion and can inherently treat problems of strong 
interactions that go beyond the pairwise summation ap- 
proximation. We use this new expression for the van der 
Waals interaction energy first to rederive the standard 
Lifshitz result for planparallel geometry as well as the 
interaction between two spheres, and then, to show its 
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versatility, we treat the problem of self-interaction of a 
uni-axial ellipsoid. 

In |7|, the nonretarded zero temperature form of the 
Van der Waals energy of a dielectric medium is given as 

F= T |^Tr (In (l + e" 1 [V, e] • VG)) = f° £/(«), 
Jo L ~ Jo L ~ 
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where e(r, r') = e(r)<5(r — r') is the local dielectric re- 
sponse function, G(r,r') = — l/47r|r — r'| is the scalar 
Green function of the Laplace operator and u is the imag- 
inary frequency that e depends on. In what follows we 
will always state out results in terms of a dimensionless 
single imaginary frequency contribution /(it) to the in- 
teraction energy defined above. The expression [, ] is a 
commutator and the trace is to be understood as setting 
the initial and final point r of the operator to be the same 
and then integrating over the space r. 

Let us furthermore assume that the space is partitioned 
into domains of constant e. For a local dielectric function 
[V,e(r, r')] = (Ve(r)) holds true and therefore 

/(u) = Tr(ln(l + e- 1 (Ve)-VG)). (2) 

The implicit function that defines the boundary S is 
given as E(r) — 0. The gradient of the dielectric function 
is thus Ve(r) = 5e 5(E(r)) VS(r) = Se D s (r) n r , where 
Se is the jump in the dielectric constant going from the 
region of negative X to positive, n r is the surface normal 
and Ds is the surface delta function such that 

J d 3 r D s (v)g(r) = j> dS g(v) (3) 

for any function g. While all the operators defined thus 
far act on the full three dimensional space, we will now 
show, that under the above assumptions they can be re- 
duced to two dimensional expressions that only operate 
on the boundaries between dielectrically homogeneous re- 
gions of space. Let us first write the trace 
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as a series, where 

T T y =D s {r) 2A r n r - VG r , r >, (5) 
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and A = |e _1 £e. While A is a function of the coordinate, 
it is by assumption a constant along any single boundary 
between two dielectrics. In a system with more than 
a single boundary we may have A's that are different 
for different boundaries. The expression for A depends 
also on e" 1 , which is however ill defined on the boundary 
itself. As will be shown below, in order to have agreement 
with the Lifshitz case, one needs to set 
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if 62 is the dielectric constant of the material on the side of 
the boundary to which the normal n r is pointing, and ei 
is the dielectric constant of the material on the opposite 
side. 

The expression for the trace of a power of T is 



Tx(T") = J dV 1 ) J dV 2 )...y 

T r (l) (2) T r (2) ,(3) ■ ■ ■ T r (n) r (l) ■ 



d 3 r (n) 



(7) 



By inserting the definition ([5]) we obtain 

Tr(T") = 2" J d 3 rM J dM 2 > ...J dM"> 

L»s(r (1) ) [A r(1) n r (i) • VG r(1 ) >rW ] x 
xD s (r (2) ) [A r (2)n r( 2) • VG^)^)] x . . . 
...xfls(rW) [V,n rW .VG r(n ,, r (i)] 
or, according to equation 

Tr(T") = 2" j dS {1) j> dS {2) . . . I dS {n) (8) 

[A r ii)n r (i)'VG r (i)j(2)] [A r ( 2 )n r (2) ■ VG r (2) ir <3)] ... 
. . . [A r („)n r( „) • VG r (») ir (i)] . 

The expressions for Tr(T n ) are thus evidently reduced to 
surface integrals. The relevant operators can be therefore 
considered to act only on the surface S and not on the 
whole three dimensional space. If we now define the main 
operator that acts between two points r, r' on the surface 
as 

K r , T > = 2A r n r • VG r , r / (9) 
defining at the same time the surface trace as 

Tr s A T>t , = IdS A r , r) (10) 



we then see that Tr T n = Tr s K n . This allows us to 
re-sum the equation J?]) as 



Tr (In (1 + T)) = Tr s (In (1 + K)) 



(11) 



and therefore formulate the interaction energy expression 
succinctly as 



f(u) = Tr s In (1 + 2A r n r • VG r 



ln|l + Ki|.(12) 



where the trace Trg now stands for an integral over the 
surface S rather than the whole space. In the last line we 
introduced the eigenvalues of the operator K, denoted by 
ki, and evaluated the trace as a sum over these eigenval- 
ues. 

To test the above approach, let us consider the Lifshitz 
case of a planar slab of thickness I composed of one di- 
electric, surrounded by different semi-infinite dielectrics 
on either side. Let us define the normal n to point out- 
wards on the surface of the slab. The operator K does 
not have any contributions arising from the interaction of 
the elements of one wall with other elements of the same 
wall, since the gradient of the Green function between 
such two elements is perpendicular to the boundary nor- 
mal. 

Let the coordinates x, y lie in the plane of the slab and 
the z coordinate perpendicular to it. Let the eigenvector 
components be denoted as Q (x, y) with j = 1,2 for the 
two boundaries. Due to the symmetry of the problem we 
may expect the solutions to take the form 



6 Q ' j \x,y) =exp(iQx)f3(j)- 



(13) 



We choose the x direction to be along the wavevector Q 
without loss of generality. Explicitly inserting the Green 
function expression into the eigenvalue equation 



(14) 



we obtain for j ^ j' 



KiP(j) = -2A, 
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dx\dyi I exp(iQa;i)/3(j') 
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Integration over y\ gives 

The integration over x\ requires contour integration de- 
pending on the sign of Q and yields the system of equa- 
tions 



Ki p(l) = -A 2 exp(-|QK)/?(2), 
Ki/?(2) = -A ie xp(-|Q|I)/3(1), 



(16) 



as we explicitly insert the two possibilities for j,f. For 
a given Q this system has two eigenvalues, namely 



(17) 



^ = i^AiAa exp(-|QK)- 
We can therefore write the /(it) (fT2|) as 

f(u) = In (1 - Ax A 2 cxp(-2\Q\£)) , (18) 
Q 

or equivalently, if Qs are assumed to fulfill periodic con- 
ditions on a flat plate of an area A, then taking this area 
to infinity, we get for the interaction energy 



F/A 



du 
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ln(l-A 1 A 2 exp(-2|g^)), 
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which is exactly the non-retarded Lifshitz result [lj|. 

In a general case, the equation (fT2")) cannot be solved 
analytically but rather requires a proper discretization 
scheme for the operator §§§ on the surface. Let us assume 
that the surface is split into a set of discrete boundary 
elements Si, along which the eigenvectors of the operator 
K are constant. Higher order schemes may of course be 
employed, but for clarity we will deal with this simplest 
example. It is worth noting that 



h Si ,r< = [ dS r n r • VG r . r , = ifi Si (r'), 
Js, 4tt 



(20) 



where fis 4 (r') is the solid angle of the surface Si as seen 
from the point r'. To obtain the matrix elements of K, 
we also average above result over the surface Sj such that 



K it j = g S% I c/.s',.' /' 



Si,r>- 



(21) 



The operators introduced thus far are not well defined 
for very short distances and lead to divergencies that are 
due to the local dielectric response assumption 8] . These 
divergencies can be regularized by a simple ansatz 
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where dij is the distance between the centers of the 
boundary elements i and j and a gives the estimate of 
the nonlocal response distance. By calculating the eigen- 
values for a discrete matrix of this type, the interaction 
energy can then be calculated according to equation (fl"2"l) 
using the approximate eigenvalues of K. 

We test this method of calculation in two cases, namely 
the interaction between two dielectric spheres and the 
self-interaction of an uniaxial ellipsoid. In both cases 
we represent the spherical surfaces via a recursive sub- 
division of an icosahedron that limits towards a sphere, 
as shown in Fig. [TJ This triangulation procedure is ex- 
tremely useful since in this case a closed expression exists 
for the solid angle Og i (r') @, which, for a triangle Tj with 
the vertices r$, reads as 

Q T . ( r ') = 2 tan" 1 [Ri ■ (R 2 x R 3 ) / (23) 

/ (R1R2R3 + R\ R2 • R3 + Ri R3 • R-i + R3 R-i • R-2)] , 

where Ri = r, — r'. Our reformulation of the van der 
Waals interaction energy was obtained exactly in terms 
of these solid angles, eq. (fT2")) . The same triangulation 
procedure can be used for any other bounding surface. 

The interaction energy calculation for a pair of spheres 
of radius R is given in Fig. O In all calculations we 
use the largest possible value of A = 1. We may see 
that the energy at large separation, that obviously corre- 
sponds to twice the self interaction of individual spheres, 
is a relatively poorly convergent function of the degree of 
discretization. The convergence can be much improved 
by simply subtracting the self interaction energy from 








FIG. 1: The triangulation representation of a sphere as used 
in calculations. We start with an icosahedron and on each 
step subdivide each face into 4 smaller triangles. This trian- 
gulation is essential and allows us to use an explicit form for 
the solid angle Qs { (r'). 




FIG. 2: The van der Waals interaction energy calculation for 
a pair of spheres of radius 1 with a = 0.2R and A = 1 as a 
function of the separation d between the spheres as compared 
to the energy of two well separated spheres. The inset shows 
the raw energy computation. The dotted lines represents the 
single, the dot-dashed line the double, the dashed line the 
triple and the full line the quadruple icosahedron subdivision 
of the two spheres. 



the total energy with the result that the difference now 
converges a lot faster. The short range kink in the inter- 
action energy is a consequence of the short range cutoff a 
that represents the nonlocal dielectric response. The van 
der Waals interaction energy of two spheres can be also 
calculated analytically via the secular determinant of the 
field modes [8j], yielding a closed form expression accu- 
rate up to 0(g?~ 12 ), where d is the distance between the 
sphere centers. In Fig. [3] we validate our result by com- 
paring it to this expression. Obviously the large distance 
behavior does not depend on the short range cutoff, as 
indeed one would expect, and agrees pleasingly with the 
analytical result Q. The short range kink is displaced 
towards progressively smaller values of the separation as 
the short range cutoff is diminished, see Fig. [31 
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FIG. 3: The van der Waals energy as a function of the sepa- 
ration for the same system as in the figure [2] but with varying 
the cutoff a = 0.2R (full), a = OAR (dashed) and a = 0.8R 
(dot-dashed). All the results are given for the third subse- 
quent icosahedron subdivision. The dotted line represents 
the analytical result ||. 
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FIG. 4: The van der Waals energy contribution for a uniaxial 
ellipsoid of constant volume as a function of the main semi- 
axis half length relative to the energy of a sphere. The inset 
shows the raw free energy computation. The dotted lines rep- 
resents the single, the dot-dashed line the double, the dashed 
line the triple and the full line the quadruple icosahedron sub- 
division of the sphere used to obtain the ellipsoid. 

The van der Waals self-interaction of a uniaxial ellip- 
soid is given in Fig. |4] to illustrate the power of our 
approach. Again, convergence of the raw value of the 



energy as a function of the density of triangulation is 
slow, but if we subtract the energy of the spherical con- 
figuration the convergence is again much faster. We see 
that the spherical configuration has the highest van der 
Waals energy and is thus unstable. This shape instabil- 
ity is driven purely by shape-dependent van der Waals 
interactions under the assumption of a negligible surface 
energy contribution. In general, the exact location of the 
point of instability would depend on the dielectric dis- 
continuity at the surface of the ellipsoid as well as on the 
surface tension and possibly on the elastic deformation 
energy [Til ] . 

In this work we presented a general numerical method 
to calculate non-retarded zero temperature van der Waals 
interactions for a class of systems that are composed of 
spatial domains of constant dielectric response. This was 
achieved by reformulating the interaction energy expres- 
sion that contains volume traces in such a way that now it 
contains only surface ones. This reformulation yields it- 
self to a straightforward numerical implementation based 
on multiple triangulation of the shapes of interacting 
bodies. We showed that the proposed method reduces 
to standard results for two planparallel semispaces and 
two spheres and that it can be used to shed light onto 
cases which have heretofore eluded an exact or even an 
approximate analysis. 

Though our approach has been formulated in the 
framework of non-retarded zero-temperature van der 
Waals interactions, it appears to us that it should be 
quite straightforward to extend it also across the com- 
plete domain of spacings between interacting bodies as 
well as to the case of finite temperatures, where the fre- 
quency integral is simply turned into a Matsubara fre- 
quency summation. 
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